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Abstract 

In this paper we investigate the solubility of a hard - sphere gas in a solvent modeled as an 
associating lattice gas (ALG). The solution phase diagram for solute at 5% is compared with the 
phase diagram of the original solute free model. Model properties are investigated through Monte 
Carlo simulations and a cluster approximation. The model solubility is computed via simulations 
and shown to exhibit a minimum as a function of temperature. The line of minimum solubility 
(TmS) coincides with the line of maximum density (TMD) for different solvent chemical potentials. 

PACS numbers: 61.20.Gy,65.20.+w 
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I. INTRODUCTION 



Solubility is the capacity of a given solute to form a homogeneous solution in a given 
solvent. The solubility of one substance in another is determined by the balance of inter- 
molecular forces between the solvent and solute, and the entropy change that accompanies 
the solvation. Factors such as temperature and pressure will alter this balance, thus changing 
the solubility. 

In the process of dissolving a solute, heat is required to break the intermolecular forces 
between the solute particles and heat is given off by forming bonds between solute and 
solvent particles. Consequently, increasing the temperature can make it easier or harder to 
dissolve solute particles in a solvent. 

In the case of dissolving solids, the energy required to break the intermolecular forces 
between the solute particles is usually higher than the energy liberated by forming bonds 
with the solvent, therefore, the ability of the solvent to dissolve the solid solute increases 
with the temperature. 

In the dissolution of gases in liquid solvents the scenario is more complex. For a low 
density gas phase, the solubility is the equilibrium constant of the process in which a solute 
particle dissolved in the liquid evaporates into the gas phase. This quantity can be express in 
terms of the inverse of the the Henry's constant. The equilibrium constant can be computed 
by the balance between the solute-solute, solvent-solvent and solute-solvent energies, namely 
W = (w so i u te- solute + w S oiv-soiv - 2w so iute- soiv) , and it is proportional to e w/kBT . Therefore, 
if W > 0, the solubility increases with temperature, otherwise it decreases. This approach, 
based in a simple "Bragg- Williams approximation", does not allow for entropic effects besides 
the combinatory term that differenciats solute from solvent particles 1 . 

This simple picture is valid for a number of solvents and solutes but not for non-polar 
gases in water. The solubility of noble gases in water decreases with the temperature until a 
certain threshold, and then it increases^"-. This experimental result has also being observed 
in the solubility of hard spheres simulated for effective^ and atomistic models of water-. 

Besides the unusual behavior of solubility, water also exhibits other thermodynamic and 
dynamic anomalies. Unlike other liquids, water does not contract upon cooling and the 
specific volume at ambient pressure starts to increase when cooled below T w 4°C-^. Ex- 
perimental data show that the diffusion constant D, increases on compression at low tem- 
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peratures T, up to a maximum, D m3jX (T) at p = PDmax{T)^ while for normal liquids the 
diffusion coefficient decreases as the system is compressed. These findings were also sup- 
ported by simulations both in atomistio^ 1 ^ and in effective models—"—. 

Even though both experiments and simulations indicate that non-polar gases exhibit an 
anomalous behavior in the solubility in the same region of pressures and temperatures where 
other anomalies are present, a clear connection between these thermodynamic and dynamic 
properties is still missing. 

In this paper we aim to shade some light into this problem by computing the solubility 
of non-interacting particles in a lattice model for water. Even though there are a number 
of two and three dimensional lattice models that would in principle exhibit the anomalies 
present in water->^£~— , here water is represented by the associating lattice gas^~— that 
have already shown the density and diffusion anomalies described above. The solute is 
represented by pure hard sphere in an attempt to model simple gases. Solubility versus 
temperature is calculated and tested for the presence of a minimum. We show that for a 
hard sphere solute, the temperatures of the maximum density and the temperatures of the 
minimum solubility coincide. This result is explained in the framework of the two length 
scales potentials and the presence of thermodynamic and dynamic anomalies in water-like 
models. 

The paper is organized as follows. In sec. II the model is presented, in sec. Ill the results 
obtained by employing Monte Carlo simulations are shown and in sec. IV the problem is 
analyzed by a Cluster approximation. Our conclusions are summarized in sec. V. 

II. THE PURE WATER MODEL - ALG 

We consider the three dimensional Associating Lattice Gas Model (ALG) introduced 
by Girardi and coworkers^. Dynamic and thermodynamic properties of the system were 
previously obtained by Monte Carlo simulations^^ 1 ^ and analytical methods^. The model 
consists of a body centered cubic lattice where sites can be either occupied (<7j = 1) by a 
water molecule or empty (c^ = 0). Besides the occupational variables there are eight arm 
variables, Tj, that describe the molecule orientation. Four arms are the usual ice bonding 
arms, with Tj = f, while the other four are inert arms, with r< = 0. The bonding and 
non-bonding arms are distributed in a tetrahedral arrangement and there are two possible 
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configurations for a water particle, as shown in Fig. [TJ 
The Hamiltonian that describes the system is given by 

H = } j OiOj (e + 7 nTj) 
(hi) 



(1) 



where = 0, 1 is the occupational variable, e is a van der Waals like energy, 7 is the bond 
energy and Tj = 0,1 corresponds to the arm variable that represents the possibility of a bond 
between two nearest-neighboring (NN) particles. A bond is formed when two neighboring 
particles have bonding arms (r = 1) pointing to each other. The parameters are chosen to 
be e > and 7 < 0, which implies an energetic penalty on neighbors that do not form a 
bond. 



Figure 1: The two possible configurations (A and B) of the water molecules in a BCC lattice. 

The ground state of the pure solvent system can be inferred from inspection of the 
equation HJ and taking into account an external chemical potential p. At zero temperature, 
the grand potential per volume is Q = e — pp, where p is water density and e = H/V. At 
very low values of the chemical potential, the lattice is empty and the system is in the gas 
phase. As the chemical potential increases, coexistence between a gas phase, with p gas = 0, 
and a low density liquid (LDL) with Pldl = 1/2 is reached, at p gas -LDL = 2(e + 7). In 
the LDL phase each molecule has only four neighboring molecules, forming four hydrogen 
bonds (HBs), and the grand potential per volume is Qldl = £ + 7 — p/2. As the chemical 
potential increases even further a competition between the chemical potential that favors 
filling up the lattice and the HB penalty that favors molecules with just four NN appears. 
At 

I^ldl-hdl = 6e + 27, the LDL phase coexists with a disordered high density liquid 
(HDL), with Phdl = 1- in the HDL phase, each molecule has eight NN occupied sites, but 
forms only four HBs. The other four non-bonded molecules contribute with positive energy, 
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which can be viewed as an effective weakening of the hydrogen bonds due to distortions of 
the electronic orbitals of the bonded molecules. The grand potential per volume is then 
&hdl — 4e + 27 — fi. Both phases are illustrated in Fig. [2] and Fig. [3j 

The non-zero temperature properties of the model were obtained from Monte Carlo grand 
canonical simulations 25 in a previous publication, for 7/e = —2. Reduced parameters were 
defined as 

rp _ k B T 

1 

e 

11 = 1 ■ (2) 

The chemical potential versus temperature phase diagram of this water model is illus- 
trated in the inset of Fig. @J Continuous transitions were investigated through analysis of the 
specific heat, of energy cumulants and of the order parameters. First-order transition points 
were located from hysteresis of the system density as a function of the chemical potential. 

At low chemical potentials and low temperatures there is a first-order transition (circles) 
between the gas and LDL phases. As the temperature is increased this transition becomes 
continuous at a bicritical point. At intermediate chemical potentials and high temperatures 
there is a continuous transition (dashed line) between the gas and the disordered phase. At 
high chemical potential and low temperatures there is a first-order phase transition (squares) 
between the LDL and the HDL phases. This transition becomes critical at a tricritical point 
(higher temperatures). This system also presents the density anomalous behavior observed 
in liquid water. At fixed pressure, the density increases as the temperature is decreased, 
reaches a maximum and decreases at lower temperatures. The line of the temperatures of 
maximum density (triangles) is illustrated in the inset of Fig. HJ 



III. NON-POLAR SOLUTE IN ALG SOLVENT - SIMULATIONS 

We investigate the solution phase diagram for inert solute particles. Non-polar solutes are 
modeled as non-interacting hard spheres. Thus each lattice site can be empty or occupied 
either by a water particle or by a non-polar solute. Simulations were run as follows: solvent 
molecules are inserted or removed via a Grand-Canonical Monte Carlo algorithm, while a 
fixed number of solute particles are allowed to diffuse throughout the lattice, following the 
Metropolis prescription. 
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Figure 2: In the ordered low density liquid phase half of the lattice is occupied by particles. Water 
(filled circles) and empty sites (open circles) 




Figure 3: In the ordered high density liquid phase the lattice is fully occupied by particles. 
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Figure 4: Chemical potential vs. temperature phase diagram for the original ALG model (inset) 
and ALG with 5% solute concentration. 

Fig. H] illustrates the effect on the water chemical potential versus temperature phase 
diagram upon the introduction of 5% solute (p so iute — 0.05 in volume fraction). At T = 0, 
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the solute will be randomly distributed in the lattice. In the gas phase, because solute is 
inert, the semi grand-potential free energy per volume is the same as the grand potential of 
the solute free case, fl gas — 0. In the LDL phase, the solute particles occupy the empty sites 
which correspond to half the lattice: for 5% of solute particles, 10% of the empty sites are 
occupied, as illustrated in Fig. Consequently, the density and the semi grand potential 
free energy per volume are the same as that of pure water, pldl = 1/2 and Qldl = £+7—1/2 
respectively. The same occurs in the phase boundary between the gas and the LDL phases 
at zero temperature, where p gas -LDL = 2. This semi grand-potential free energy per volume 
is only correct if the solute occupies up to 50% of the lattice. 

In the HDL phase the addition of the solute changes the density and the grand potential 
per volume when compared with the free solute case. At zero temperature, in the solute free 
case, the lattice is full. As solute is added the particles replace solvent particles, changing the 
energy, as can be seen in Fig. [6j The solute can occupy a number of different configurations 
but it will choose the structure with the lowest grand potential free energy. In order to 
identify which is this configuration we compare two cases: all the solute forming a large 
cluster and the solute distributed randomly but not occupying adjacentes sites. In both cases 
Phdl — ^-—p solute but the grand potential free energy is given respectively by ^5/D^ r = 2(3e+ 
2<y)-(Ae + 2<y)p so i ut e- p{l- psoiute) and by VL r ^ d £ m = 4e + 2>y- (Qe + 2>y) p so i ute - p(l -p so i ut e)- 
The last case is the scenario with lower grand potential free energy, therefore Even thought 
the density and the grand potential per volume is different from the free solute case, the 
phase boundary between the LDL and the HDL phases is given by the same value as that 
of pure water, namely Pldl-hdl = 6e + 2j. 

At low temperature, and T^O, the phase boundary between the LDL and the gas phase 
depends on the entropy related to the addition of the solute. In the gas phase the solute can 
be located at any site, while in the LDL phase the solute goes only to the empty half-lattice 
(see Fig. [5]), which makes the entropy of the gas phase larger than the entropy of the LDL 
phase. This stabilizes the gas phase, moving up the chemical potential of the gas-LDL phase 
boundary with respect to the free solute case. 

As to the LDL-HDL phase boundary at low T ^ 0, distorion of bonds in either phase is 
small, so that at very low temperatures no entropic difference between the two phases due 
to the presence of the solute is observed. However, as the temperature is increased and the 
bond network becomes more disordered, the entropy increment caused by solute is larger in 
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the LDL phase, thus shifting the phase boundary to higher chemical potentials. 

The critical lines illustrated in Fig. H] were obtained from the peak in the specific heat 
and through Binder cumulant methods^^. The critical line between the disordered phase 
and the HDL phase, at high temperatures (r line), and the critical line between the HDL 
phase and the LDL phase (A line) are similar to the lines obtained in the solute free case^. 
The shift in these lines to lower temperatures is due to the higher entropy of the solute in 
the disordered phase. 

Besides the A's and the r critical lines present even if no solute is present, the system 
exhibits a new critical line, which we call 77, at high chemical potentials. In order to inspect 
the new phase, we investigate the solute spatial distribution through its radial distribution 
function, g(r), given by 

(Ei=i Ej>i o-i<W, [Ti - r^Doo 

where the mean values denoted by < ... >^ are obtained at constant temperature for T = T 
and T — >■ oo, and 5 is the Kronecker's delta. 

Figure [7] illustrates the radial distribution function for 5% of solute for two different 
temperatures, T = 0.40 and T = 0.70. At high temperatures the radial distribution is 
characteristic of a disordered, fluid-like system, while at low temperatures solute presents 
an ordered (solid-like). 




Figure 5: In the ordered low density liquid phase half of the lattice is occupied by particles and the 
other half contains solute. 

To further analyse the structure of the new phase, we have measured sublattice solute and 
water densities. The lattice is divided into eight sub-lattices, as shown in Fig. [HJ The result 
is shown in Fig. [9] for chemical potential fixed at p, — 3.00. At high temperatures, both 
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Figure 6: In the ordered high density liquid phase the lattice is completely occupied either by solute 
or water particles. Solutes are the bigger circles. 
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Figure 7: Solute - solute radial distribution function, g(r), for solute concentration 5%. From top 
to bottom T = 0.40 and T = 0.70, chemical potential p. = 3.00. At low temperature g(r) exhibits 
some structural order while at high temperature g(r) is liquid like. 

solvent and solute are randomly distributed, with no preferentially occupied sub-lattice. As 
temperature decreases and the vertical critical line is crossed, the solute particles become 
organized in four sub-lattices, while the solvent fills the remaining space in the eight sublat- 
tices. Therefore, the 77-line separates the higher temperature homogeneous fluid phase from 
a dense phase in which the solute organizes itself, named SLDL/HDL phase. It is not clear, 
however, if inside each sub-lattice the solute particles are randomly distributed or if they 
cluster, forming a high density solute and a low density solute region. The behavior of the 
radial distribution function suggests that the later is the case. This issue is checked through 
the cluster approximation of the next section. 
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Besides the critical lines and the coexistence lines, this system with and without solute, 
exhibits an anomalous behavior in the density. As the temperature is decreased at constant 
pressure the density increases, and at very low temperatures, it decreases, thus presenting 
a maximum. The line of maximum density temperatures at different chemical potentials is 
illustrated as the line with triangles in Fig. HI Due to the entropic effects the TMD line is 
located at higher temperatures when compared with the solute free case. 




Figure 8: The eight sub-lattices division. 




Figure 9: Density of sub-lattice as a function of temperature for chemical potential /i = 3.00 and 
solute concentration 5%. Left graph shows solvent molecules and the right one, the solute. 

Having explored the effect of introducing inert solute upon the solvent chemical potential 
versus temperature phase diagram, we investigate the solubility of these non-interacting par- 
ticles. The solubility is derived from the solute excess chemical potential which is computed 
through Widom's insertion method namely 

fj, ex = -T \n(e-^ Me ) Q . ( 4 ) 
fO 



where e so iute represents the energy increment due to the added solute. Note that the solutes 
interact only via excluded volume. In the limit of low solute concentration, the solubility is 
then given by 



Figure [10] illustrates the solubility parameter T as a function of temperature for solute 
concentrations 5% for various chemical potentials. It can be seen that the solubility exhibits a 
minimum for a certain range of chemical potentials. The temperatures of minimum solubility 
(TmS) are shown with star symbols. The TMD and the TmS lines coincide. This coincidence 
can be understood as follows. 

The solubility is related with the amount of energy required to include a solute particle in 
the system. In the gas phase, as the temperature is increased, the solubility decreases because 
the density of the gas phase increases with temperature making more difficult to include an 
extra solute particle. In the HDL phase the density of solvent decreases with the increase 
of temperature, therefore the solubility decreases with the decrease of temperature. In the 
LDL phase a very interesting behavior is observed. The density, which at zero temperature 
is 1/2, increases with the temperature. In this low temperature region the increase of 
the solvent density leaves less space for including a new solute particle and the solubility 
decreases. However, as the temperature is increased even further, the density of the LDL 
phase reaches a maximum and decreases due to entropic effects. Consequently, the solubility 
reaches a minimum at the maximum of solvent density, and increases with T as the system 
becomes less dense. The decrease of the solubility with the increase of the temperature is 
also observed in a minimal lattice models proposed by Widom and collaborators^. The 
focus of these different studies on the model properties have been on the relation between 
strength and range of the solvent mediated attraction, thus density effects, which would 
promote increasing solubility at larger temperatures, were left aside^ 1 ^. 

Obviously the solubility governed solely by density effects is only possible because the 
solute particles, in our study, are non interacting. 




(5) 
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Figure 10: Solute excess chemical potential vs. temperature for different values of the water chemical 
potential. 

IV. NON-POLAR SOLUTE IN ALG SOLVENT - CLUSTER VARIATIONAL AP- 
PROXIMATION 

In this section, we present results obtained from the Cluster Variational Approximation, 
introduced for the pure solvent version of this model by Buzano et al.—. The approximation 
is based on the assumption that an elementary tetrahedron cluster of sites, connecting the 
four sub-lattices (Fig. EJ) is thermodynamically representative. 

The energy per site is written as 

W = PijklHijkl ■, 

i,j,k,l 

where Pijja is the joint probability for the configurations of the cluster in the four sub-lattices 
1, 2, 3 and 4, and Hijhi represents the Hamiltonian for the cluster. 

For interactions restricted to nearest-neighbors the Hamiltonian Tiijki can be decomposed 
in the form 

Hijkl = TXij + Hjk + H-ki + Hu > 
where the pair Hamiltonian can be written as 

Hij = eaiaj - -o-iCTjh,: - — — . 

Here, G\ = 1 if site i is occupied by a water molecule and U{ = 0, otherwise. The new 
state variable b stands for solute occupation, with 6j = 1 if site i is occupied by a solute 
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particle and b{ = otherwise, hij represents bond states and is equal to 1 if there is a bond 
between the two waters at sites i and j, and equal to zero for non-bonding waters. To the 
parameters of the previous section, hydrogen-bond interaction 7, non-bonding penalty e > 
and solvent chemical potential fi we add the solute chemical potential v. 
By minimizing the grand-canonical free energy 



where pj = J2 jkl p ijk h v) = £ iW Pyw, v\ = Y,ijiPw and Pi = ^, ; /, /'<,/.•/• According to the 
natural iteration method^, for fixed values of T, e, 7, fi, and u, one may different quantities 
of interest, such as solvent and solute concentrations, hydrogen bond fraction or sub-lattice 



In order to compare the cluster approximation results to previous simulation results, 
we search for the desired solute concentration (p so iute) by varying the corresponding solute 
chemical potential v. In certain ranges of the thermodynamic parameters, hysteresis loops 
were found. 

The phase diagram for 5% solute concentration {p so iute = 0.05) is displayed in Fig. [TTJ 
It is qualitatively similar to that of simulations, and the new phase SLDL/HDL in which 
solute is ordered is present. In that phase, a solute droplet occupying a sub-lattice forces 
the neighboring solvent to order as a LDL. All other non solvent molecules are in the HDL 
phase. In other words, as the temperature is lowered and one crosses the fluid - SLDL/HDL 
line for a given value of /1, a second-order phase transition takes place, and the solute, which 
was dispersed in the fluid phase, begins the formation of sub-lattice clusters, intercalated 
with LDL water. Thus, the SLDL/HDL phase is a region of coexistence between clustered 
and dispersed solute, or between LDL+solute and HDL-no solute. 

The picture we propose for the SLDL/HDL phase is supported by behavior shown in Fig. 
[T2| for solute density as a function of solute chemical potential z/, at fixed T and solvent 
chemical potential /z, in the vicinity of the HDL - SLDL/HDL line. As v increases, the 
solute density goes through a discontinuity. The diagrams show coexistence between a low 
density, with solute dispersed in HDL solvent, and a high density, with solute grouped in 
LDL solvent. Maxwell construction for the curves of Fig. [T21 yield the densities of the two 
coexisting phases shown in Table 1. It can be seen that for T = 0.4, for both values of 
the solvent chemical potential, as well as for T = 0.7 and fi = 2.4, the case of 5% solute 




densities. 
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Figure 11: Solvent reduced chemical potential p versus reduced temperature phase diagram. The 
Gas-LDL and LDL - SLDL/HDL coexistence phases, critical and TMD lines are shown. 




V 



Figure 12: Solute concentration p so i u t e as a function of the chemical potential v for some values of 
temperature and water chemical potential p near the fluid - SLDL/HDL transition line. Dashed 
lines give the coexistence lines. 

density (p so iute = 0.05) is between the low and high density limits (pi and ph). The three 
cases correspond to points on the left of the fluid - SLDL/HDL 77-line, although very near 
to the line in the case of T = 0.7. On the other hand, for \x = 3.0 and T = 0.7 (clearly on 
the right of the HDL - SLDL/HDL 77-line), the low density at coexistence is p iow = 0.078, 
larger than 0.05 implying solute must be homogeneously distributed in the solution. 
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Table I: Solute chemical potential at the coexistence u coex and the respective concentrations of 
solute in both phases. 

V. CONCLUSIONS 

In this paper we have analyzed the effect of adding non interacting solute particles to the 
Associating Lattice Gas Model for water. We have shown that the solute particles change 
the chemical potential versus temperature phase diagram substantially. In the first place, by 
enhancing entropic effects of the disordered and gas phases, addition of solute displaces the 
coexistence and critical lines accordingly. Moreover, a new phase appears, with coexistence 
between HDL and LDL water induced by the presence of solute. Thus, in addition to the 
A and r critical lines present in the solute free case, a third order-disorder continuous 77 
transition at which the solute orders itself is present. 

The solubility of this non interacting system was also computed and we have found that, 
as the temperature is decreased for fixed solvent chemical potential, the model solution 
exhibits a minimum in the solubility. The temperature of minimum of solubility coincides 
with the temperature of maximum density. This coincidence is due to the fact that the 
solute is purely hardcore. For hydrophobic or hydrophilic solutes this coincidence should to 
be lost. 

Our results suggest a very simple mechanism for understanding the minimum of solubility 
in simple gases. 
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